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ABSTRACT 



Recent interest in military VTOL/STOL aircraft employing unprepared 
landing sites has led to interest in the problem of landing surface ero- 
sion. Surface erosion is caused by the aerodynamic forces on ground 
particles existing within the flow field of an impinging jet. The 
inviscid flow field is discussed and the viscous ground boundary layer 
is analyzed utilizing both theory and available experimental data. A 
mathematical model of the process of entrainment of ground particles is 
constructed. Erosion rates in the form of erosion profiles are pre- 
dicted for selected jet configuration and types of terrain. A criterion 
for entrainment, due to both lift and drag, was found and presented for 
selected distances from the jet centerline. 
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NOTATION 



O 

erosion rate - ft/sec 



coefficient of drag 
coefficient of skin friction 
nozzle diameter - ft, 
particle diameter - in. 
force - lb. 

lift force due to the effect of the ground - lb. 
lift force due to shear - lb. 
mass - lb sec^/ft. 



shear velocity - ft/sec. 

radial distance from jet centerline - ft. 

radius of particle - in. 

kinetic energy - ft lb. 

velocity parallel to the ground - ft/ sec. 
velocity perpendicular to the ground - ft/ sec. 
coordinate parallel to the ground 
coordinate perpendicular to the ground 

vertical distance at which ^ ^ UmcKK in turbulent boundary layer 
height of nozzle above the ground - ft. 

r\ A 

density - lb sec^/ft . 
viscosity - lb sec/ft^. 



rotational flow parameter = 





number of strips used in strip analysis 
static pressure on the ground Ib/ft^. 



© 



angle measured from center of particle 
vor ticity 



COo 

S boundary layer thickness 

Subscripts 

crit. value at boundary layer transition 

J indicates value associated with jet 

L value at local point on the particle 

a property of air 

I value in inviscid flow 

p value associated with a particle 

c measurements made on cylinder in strip analysis 

s static 

A ambient 

n measurement made on a strip 

w value associated with effect of wall constraint 

S value associated with shear effect 
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1 . 



Introduction. 



The increasing interest in, and use of, VTOL aircraft and other 
vehicles with vertically directed fans and jets, particularly in mili- 
tary operations, has brought about concern for the deterioration of 
unprepared or soft landing sites. The downwash impingement caused by 
these aircraft results in a radially developed viscous flow field on the 
ground, producing forces adequate to cause entrainment of the ground 
particles. The reasons for concern with this problem are numerous. In 
addition to possible damage to the aircraft itself and the erosion of 
the landing surface there are the accompanying hazards to ground per- 
sonnel, damage to ground support equipment, poor visibility afforded the 
pilot, and the violation of the military concept of concealment. The 
present study was designed to construct a mathematical model of the 
ground erosion process to allow prediction of the initial erosion rates 
for various types of terrain for various nozzle velocities and geome- 
tries. Previous work Q, 14, 15, 22, 23, 25]] on the problem of erosion 
due to downwash impingement has been primarily an experimental examina- 
tion of the problem with small scale jets. 

The critical mechanism in the erosion process is the interaction of 
the ground boundary layer and the terrain particles. An analysis of 
this mechanism requires an understanding of the viscous mixing problem 
of a finite free jet, with the inviscid flow field extending radially 
from the jet center line as a result of the presence of the ground plane. 
An analysis of the phenomenon which produces entrainment of the ground 
particles is basically a study of the aerodynamic forces which exist in 
the ground boundary layer, as they are affected both by the decaying in- 
viscid flow field above the boundary layer and by the constraint of the 
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ground surface. In analyzing the boundary layer forces available theory 
was utilized whenever possible, but since the radially developed bound- 
ary layer runs the gamut of stagnation point flow to turbulent flow and 
then complete decay, there are areas for which no theory or exact solu- 
tion exists. To fill these holes, experimental data were used. Where 
exact solutions were used for the velocity profiles within the boundary 
layer, verification was made with experimental results if data existed. 

The analysis was restricted to single jets of uniform dynamic pres» 
sure loading and circular cross-section under "no wind" conditions shown 
diagrammatically in Fig. 1. An axi-symmetric boundary layer was assumed, 
acting over a terrain of uniform particle size and density. Table I in- 
dicates the various combinations of nozzle characteristics and radial 
stations that were analyzed. The prediction of erosion rates is neces- 
sarily restricted to incipient motion, since the secondary effect of 
random collision is not readily treated by analysis and because the con- 
stantly changing shape of the ground surface is not accounted for as a 
continuing process. Although this work is restricted to consideration 
of uniform jets only, previous experiments by Morse t.23j indicate that 
the dynamic pressure distribution due to downwash from rotor blades and 
ducted fans is similar to that for a uniform jet. The uniformity of 
the jet or downwash appears to be a critical factor only for ratios of 
jet altitude to diameter, Z/D^, less than unity. 
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2. History of the Problem. 

The classic problem of an impinging jet on a normal surface is not 
new, but its application to the field of VTOL aircraft, of either rotor 
or jet type, has been generated only within the past several years. The 
problem of surface erosion resulting from impinging downwash first re- 
ceived attention from Kuhn D.4]. Kuhn’s work was experimental employing 
small scale test equipment to make dynamic pressure distribution surveys. 
An experiment was conducted by fixing the nozzle dynamic pressure and 
raising or lowering it over terrain of known condition. Testing was 
done with various types of sod, dry and wet dirt, dust, and sand. Simi- 
lar experimental tests have been conducted more recently by Morse \ 23 ~]. 

The most recent work done on this problem is that of Brady and Lud- 
wig [26]. This work consists of extensive experiments to investigate 
the characteristics of the impinging uniform jet, including both the 
inviscid flow and the ground boundary layer. Measurements of boundary 
layer velocity profiles were compared with theory and showed agreement 
to a limited degree. Some of these experimental results have been used 
in support of the present work. Vidal \ 28 ] has given an insight into 
the aerodynamic forces existing within the boundary layer. 
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The Flow Field. 



In order to consider the mechanism of particle entrainment and ground 
erosion, a thorough knowledge of the flow field of the impinging jet is 
required. This flow field can effectively be separated into three re- 
gions. The first of these regions exists regardless of the presence of 
the ground surface, and is the region of viscous decay immediately exter- 
nal to the jet nozzle. Viscous decay of a free jet is a classical prob- 
lem under uniform conditions. This region is characterized by the flar- 
ing out of the jet as the jet boundary mixes viscously with the stationary 
field external to it. Introducing a ground surface into the flow field, 
normal to the jet, does not invalidate the quantitative solution of free 
jet decay, but restricts its extent of usefulness. Experimentation has 
shown that the theory of viscous decay of a free jet is adequate for 
the impinging jet problem for nozzle heights of greater than about eight 
jet diameters. Since in this work, the interest is in nozzle heights 
of S , an additional region of viscous mixing with the station- 

ary boundary must be studied and determined, beginning with the point at 
which the presence of the ground surface has altered the free jet charac- 
teristics. Further study must then be made of the viscous decay of the 
jet after impingement on the ground plane. To date no satisfactory 
theoretical method exists for dealing with the entire flow region in 
three dimensions. This problem may be partially overcome, however, by 
including the viscous decay region in describing the inviscid flow field. 

Knowledge of the inviscid flow field is necessitated by the require- 
ment that the inviscid velocity existing at the upper edge of the ground 
boundary layer be known. The difficulty in finding a theoretical solu- 
tion for this region lies in the fact that although the boundary conditions 
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are well defined, the location of the free streamline is not known. Three 
dimensional theoretical solutions of this problem are limited to approxi- 
mate methods. One of these, by LeClerc ^lO] , uses an electrolytic analog 
to determine the shape of the free streamline boundary. In this method 
the solution of the inviscid flow region was accomplished by relaxation 
techniques, predicting velocity and pressure distribution. Such a 
method of solution could be made to approach the exact solution to any 
desired degree of accuracy. 

In the present work an exact solution of the inviscid flow field 
was not required. It was assumed adequate to use existing experimental 
data which compared favorably with the limited theoretical analyses avail- 
able. In most instances it was found that the theoretical solutions were 
not in good agreement with the experimental data for The approx- 

imate methods of predicting the characteristics of the inviscid flow re- 
gion by use of an equivalent inviscid jet of greater diameter and decreas- 
ing dynamic pressure, or by replacing the jet with a cylindrical vortex 
sheet of constant radius extending to the ground, have proven unsatis- 
factory when compared with experimental data. In view of this, available 
experimental data were used in this work as a solution to the inviscid 
flow field. Fig. 2 shows the result of empirical equations fitted to 
these data. 

The third flow region is that of the ground boundary layer. A thor- 
ough knowledge of this region was most important to this work as the ground 
particles are predominately immersed in the boundary layer, and the mechan- 
ics of particle movement originate with the aerodynamic forces resulting 
from the characteristics of the boundary layer. 
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4. The Ground Boundary Layer*, 

The key to the entire analysis of ground erosion is the mathematical 
model used for the boundary layer. Sound theoretical analysis of the 
boundary layer was used in this study when possible. Experimental 
results were used whenever the theory was non-existent or was in large 
disagreement with these results. For effective analysis the boundary 
layer may be divided into four regions. The first is the stagnation 
area boundary layer existing adjacent to the jet center line and extend- 
ing to .5 or to a station directly under the original free jet 

boundary. This area is laminar for the conditions investigated. The 
solution of this region utilized the Himenez solution Lisl for the uniform 
impinging jet on a perpendicular wall. Experimental results were not 
available for this region, so that the validity of this solution is not 
verified but was assumed. The Himenez solution in three dimensional 
flow, developed by Horaann [l8]j is an exact solution utilizing the Navier- 
Stokes equations for axisyrametric flows. The solution was developed 
explicityly for stagnation conditions and the velocity profiles result- 
ing from this have been employed to 0.5. The nondimens ional 

velocity profile is shown in Fig. 3. 

The second region is the laminar boundary layer extending from the 
edge of the original jet boundary to the radial station at which transi- 
tion begins. A complete method of solution for this laminar region, by 
Smith is well verified by experiment for ^/ok 4 - *5 1263, but 

gives little insight as to the limit of the laminar region and the be- 
ginning of transition. This point has been approximated by Brady and 
Ludwig £26] from considerations of neutral stability of the laminar bound- 
ary layer. A Reynolds Number at which the boundary layer becomes neutrally 
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stable can be determined [,18^ and may be converted to a jet nozzle 
Reynolds Number as a function of C|^/Dps,\^ CRiX * this basis the 

boundary layer becomes unstable at 1*0 for the jet velocities 

and nozzle diameters analyzed in this work* In view of this, and the 
fact that the invisced velocity generally reaches a peak in the vicinity 
of (Fig. 2), the beginning of transition was taken to be 

= 1 . 0 . 



The transition region is even more difficult to define. Since 
fully turbulent flow exists theoretically when the pressure gradient on 
the ground surface is essentially zero, it was assumed that the transi- 
tion region extended to the station where the pressure gradient could be 
taken to be negligible. Fig. 4 shows the static pressure distribution 
near the ground taken from an electrolytic analogy of the entire flow 
field The station at which the gradient, / (3 X > 

was taken to be zero was = 2.0. For the laminar and transition 

regions of the boundary layer experimental results were used to form the 
velocity profiles. This data, collected by Brady and Ludwig [26^ for a 
limited range of jet velocities indicate that the non-dimensional 
velocity profile was almost completely independent of the jet velocity. 
The measurements were taken at nozzle heights in the range 
to 4.0 at four locations and for the radial range of to 1.33 

at six different stations. These twenty-four configurations are assumed 
to be fairly representative of the laminar and transition regions and 
are shown in Figs. 5a through 5h. For purposes of computation, curves 
were fitted to each group of data. 

The fourth region in the boundary layer is the turbulent regime. 

It was assumed that this region begins at 2 as this is the 
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station at which the pressure gradient is assumed negligible, A limited 
amount of work has been done on theoretical solutions of turbulent bound- 
ary layers. One of these that appears applicable to the impinging jet 
is that of Glauert This single solution for velocity profiles in 

the turbulent boundary layer, assumed valid to 4. which is the 

radial limit of the analysis, has been verified by experimental results 
at 1.33 [26], The uncertainty as to the characteristics of flow 

decay limit further extension of the analysis, 

Glauert set up the boundary layer equations and found a similarity 
solution in which the form of the velocity distribution across the jet 
does not vary along its length. The solution is dependent upon the 
assumption of an effective eddy viscosity as required to satisfy the law 
due to Blasius for flow in a pipe as well as Prandtl*s hypothesis for 
free turbulent flow. The first solution is considered to be valid near 
the ground surface, and the second to be valid above some determined dis- 
tance from the ground surface. Fig. 6 shows these results in the form 
of a velocity profile, for the radial distances and nozzle heights con- 
sidered in the present work. It was found that the shape of the velocity 
profile was essentially independent of both the radial distance and the 
nozzle height, but is dependent upon the nozzle Reynolds Number. 

The turbulent region is characterized by a velocity profile which 
approaches a maximum, U r^o^x > with increasing y and then decreases with 
further increase in y. As is shown in Fig. 1 it is expected that the 
turbulent boundary layer will continue to grow as long as there exists 
an inviscid upper boundary. Since very little attention has been given 
to this phenomenon of flow decay, and because of the nature of the tur- 
bulent solution, a fictitious boundary layer thickness is assumed for 
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this region. This is defined to be the thickness that exists at 
U = Urr\(xx * found that the nature of this assumption does not 

affect the overall results appreciably. 
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5. Analysis of Aerodynamic Forces Within the Boundary Layer. 

The ground boundary layer is a non-uniform flow that is character- 
ized by a large velocity gradient or shear layer extending over a large 
portion of the thickness. In addition, the ground surface or wall bound- 
ary at the base of the flow forms a constraint on the streamlines as 
they attempt to pass beneath an obstacle bedded in the ground surface. 

There are three distinct forces that exist within the boundary 
layer. These are: 1) Drag due to the horizontal velocity component of 

the flow; 2) lift due to the proximity of the wall; and 3) lift due to 
the velocity gradient in the flow. There is no theory available which 
deals with this problem in general. Therefore it was necessary to piece 
together an approximate theory for each of the various forces. 

Drag Force. 

The drag force was the easiest to account for. To simplify the 
analysis in a non-uniform flow, a strip analysis was made in which the 
spherical particle was divided into a number of layers and it was assumed 
that a uniform flow acted upon each layer. Knowing the velocity of the 
flow on each layer, from the velocity profile, the force on each layer 
and the summation of horizontal forces were developed. The calculation 
of the total drag force employed the drag coefficient of cylindrical 
bodies measured by Hoerner 



Nif^ < i oo 
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Np; > 10^ 



Co = I, l8 



In addition to the drag force the skin friction force was calcu- 
lated for each layer and summed for the particle* The friction force, a 
result of a minute boundary layer effect on the particle itself, was 
relatively small but in some cases significant. For spheres in low Rey- 
nolds Number flows the skin friction coefficient can be approximated by: 



The tendency of the velocity gradient to cause a moment on the particle 
was neglected. 

Lift Due to Wall Proximity. 

Theories accounting for the effect of a wall constraint in non- 
uniform flow do not exist, but this effect was estimated using the exist- 
ing theory for uniform flow The effect of the ground surface con- 

straint is to distort the flow over the sphere, crowding the streamlines, 
causing a negative lift force. The theory is based on kinetic energy 
considerations, maintaining equilibrium on a sphere in the presence of a 
wall. For a sphere, in a flow parallel to the wall, an upward force is 
required to maintain equilibrium, indicating that the sphere must be 
attracted to the wall. 

The method of calculation of this lift force was similar to that of 
drag. The strip approximation method was used to determine an average 
local velocity over a particular strip. It was then assumed that this 
average local velocity was that existing in the vicinity of the stagnation 
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point, realizing that the stagnation point is shifted toward the region 
of higher velocities. The lift then becomes a function of the sphere 
size and the average local velocity, for a particle bedded on the ground 
plane. The kinetic energy of a sphere in a moving fluid near an in- 
finite fixed wall can be approximated by a first order solution £l2^: 

T- '/K C ^ ) 

where : 

A = ((ii. rr,o_ C I <■ Vyi ) 

S' <-Vyi) 

mass of fluid displaced by the sphere. 

From Lagrange's equation for equilibrium: 

^Vsx. 

If the sphere is to be maintained in equilibrium, the result is a force 
exerted on the sphere: 

- ‘'/d-t C^^dx) - 
Fy = “'/y-t C^’Vay ) - 

Making the substitutions, the force exerted on the sphere in the y or 
vertical direction becomes: 

Fy = ^ Viy 
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= AVt, - C- 'AVl."' V/g Vy ^ 



— Vs;?. V^ '' ) 

Assuming that V, V <A<^ |. 

Fy = 

where : 

^c^^VaTTri 

we get: 

p/ = ®//6 ^o. rr <‘Vyi- 

Assuming that the particle is very close to the ground such that; 

y - r << 1 
or 

y ^ r 

the lift due to the wall constraint becomes: 

rrr^ 

For the strip analysis, summation yields: 

?a.rrr^l ^-r^)/n 

Lift Due to Velocity Gradient. 

There have been a number of theoretical studies of the effects of 
velocity gradients or shear on aerodynamic forces. Most of these, how- 
ever, deal either with an unbounded flow or two dimensional flow. In 
comparing two dimensional solutions with three dimensional solutions it 
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was noted that the difference is a term in the three dimensional solu- 
tion describing the lateral component of flow. For the axi-symmetrical 
case dealt with here, this component was assumed negligible. A solution 



ployed in conjunction with a strip analysis. 

In the strip analysis, the particle was divided into layers of cir- 
cular cylinders situated such that their centers formed a line perpendicu- 
lar to the flow. Beginning with the final result of Murray and Mitchell, 
the dimensionless shear velocity on each cylinder was found to be: 



for non-uniform flow on cylinders by Murray and Mitchell 













where : 



CO - '^/c/(K L CoO ] I 










r-0 r" i C’^ 



{_ lo^ L I) ^ t Cn-tr^O]} 
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and 



CrO^"^' 







'f' Ql']) = Euler's Constant 
= 0.57721567 



The dimensional shear velocity is given by: 

‘b- 

which is of the form (when squared) : 

+C>^ ^ Z/^Q> t * 2^ Q>C.) 

Where : 



F\^ - © [ cosV-^’' c sin © - V/vi <^os\n c^sm©) 



SmV-^ C^Sm©^ VfNj^ SinV^^C. sin ^)1 

5-1, 



= 4 [^2 CriT" 



^Tvn^i CO 



X n -f 1 CO] [^s > n C.^ I ) ©J 



n-o 



.2 AB ^ Vn ^ 'T O') Qo s 2. n 0 ] • 

s i© © QcosK c •i m ©) - '/f^ s i*^K O^ ' 

ZAC. - 4[S^Cr iT" K';^^,|C0 Sin^ 2 .A^^e j • 

Sin© [^CosK C5in s')- '/rV S>nV^ C. S l n ©^"] 
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Uu is the local velocity at the center of each cylinder, and I, K and 



of the second kind respectively. N is a flow parameter dependent upon 
local velocity, radius of the cylinder and the vorticity at the stagna- 



the flow parameter was approximated since the stagnation point on each 
cylinder could not be known in advance of the calculations. In order to 
keep the shear velocity in the direction it was known to flow, a lower 
bound of N “ I could be established. If and U are evaluated at 
the geometric centers of the cylinders, experimentation has indicated an 
upper bound of N After analyzing the computed results of 

shear velocity with values of N ranging from 1.1 to 10, a value of 
N = 1.5 was accepted as being representative of the physical model. 

Having found the shear velocity the lifting force due to shear 
effect was calculated for each cylinder and summed over the particle. 

The resulting expression is: 



K* are Bessel functions of the first and second kind, and the derivative 



tion point, and is defined 




. The determination of 
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where : 



radius of cylindrical strips 
b ■= width of strip 

n = number of strips 

Noting that; 

XTT 



1 



We have: 

^ /I + 2F\<:^ +;?. c > q ) ' 

o 

s\'f^ e 

In evaluating the integral we have; 



K';,^ col 



J [f J ne s.noc'&l. 

r / Siv-\ c;^2.r\+V) © siv^er/el •= O 

■- n-o V, i 



Since the last two Integrals are identically zero. Further we may write 



■ Ztt 



f^" 

j A^Si'n^f/©^-*- { Si>^'^©-cosK^CS)-n©)f/© 

'o 'O 

~ ^/N f ■S ' © Cos K C. s 1 o S I o\^ ' n o') r / © 

•a 

^Z-TT 

'*' /n^ ■S’lio^© sioV>y‘^ C^stn ©"^cV© 
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P\^S\r\ec^Q = - Vh j -S I & cosV-\ CSn^ ' 



Sn^v^ S>n©) c/& 



Since the first and last integrals are identically zero. 

Due to the complexity of the above expressions and the difficulty 
of integration, it was desired to check this method with another method 



The alternate method of calculating lift assumes that the previously 
discussed lift due to the ground surface constraint is negligible and 
that the lift forces on the particle can be estimated by using the solu- 
tion of Hall 8 for a sphere in a two dimensional uniform shear flow. 

In addition, the approximate method assumes that the sphere rests on a 
bed of spheres of the same diameter and that the pressure on the bed is 
the ambient stream static pressureo The lift coefficient then becomes 
simply a function of the local slope of the velocity profile, , 

the local velocity on the sphere, and the radius of the sphere. In view 
of the many restrictions placed on this solution it was used only as an 
order of magnitude check on the other solutions. Fig. 7 shows a compari- 
son of the two solutions. 




of a much more approximate nature for calculating the lift due to 



shear. 
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6. Particle Entrainment and Incipient Erosion Rates. 

Only a few of the possible types of terrain to which a mathematical 
analysis of the erosion rates might be practically applied, were investi- 
gated. Uniform particles ranging in size from .003 in. to .125 in. in 
diameter were investigated. The smallest sizes represent loose dirt 
and the larger diameters represent sand or sandy gravel, with .125 in. 
being representative of small gravel type rock. It was assumed that all 
terrain was without any moisture content. The densities of the particles 
ranged from 75 Ib/ft^ for loose dirt to 125 Ib/ft^ for gravel. Table II 
lists the various terrain particles analyzed. 

There are three possibilities for entrainment of the ground particle. 
The first is that it may have enough lift force on it to be lifted directly 
from the ground surface into the free stream. The second is that after 
rolling motion is started due to the drag force, it may be able to gain 
the additional lift required for entrainment in the free stream. The 
third is that with many particles rolling at different velocities the 
random collision of particles that follows will bounce some of them into 
the inviscid stream. Only incipient motion is within the scope of this 
work, hence it was assumed that any particle which was moved contributed 
to the erosion rate. 

The net drag force and the net lift force, including both mechanisms 
of lift, were computed for each of the particles, at each of thirteen 
radial stations on the ground surface for each of the nozzle configura- 
tions studied. To calculate the initial net drag force, the previously 
mentioned ground surface model of nested particles was assumed using a 
coefficient of static friction of 
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7. 



Results and Conclusions. 



The particle size criteria for lifting entrainment and drag entrain- 
ment are shown in Figs. 8, 9, and 10 in which the particle size is plotted 
versus a load parameter at various stations with jet velocity and jet 
height as parameters* The load parameter is defined as the particle den- 
sity divided by the dynamic pressure in the free stream. The lift en- 
trainment criterion predicts that an optimum particle size, the particle 
size most prone to entrainment, exists for all radial stations and that 
the critical density increases with radial distance to a certain . 

The optimum particle size for lift entrainment occurs because particles 
smaller than the optimum are affected by much smaller velocities, and 
particles larger than optimum are affected more by a flow region where 



The drag criterion indicates that extremely small particles of 0.02 
in. diameter or less, up to very high densities will roll on the ground 
(see Fig. 8). As the particle size is increased the possibility of 
rolling decreases rapidly up to a particle diameter of 0.03 in. at which 
size the tendency to roll increases again, since the larger velocities 
in the boundary layer begin to act on the particle, until an optimum 
diameter of approximately 0.06 in. is reached. For radial distances 
less than ~ •!> velocities in the boundary layer are not suffi- 

cient to move large particles of high density. The range of optimum size 
particle for drag entrainment of approximately .03 to .09 is consistent 
for all radial stations. It is of particular significance that the range 
of optimum particle diameter for both lift and drag is very nearly equal 
to the boundary layer thickness at the particular radial station being 
considered. 



the shear gradient. 




, approaches zero. 
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The incipient erosion rates for selected densities, particle sizes, 
and nozzle configurations are shown in Figs. 11 and 12. These are found 
by combining as vectors the incipient horizontal and vertical accelera- 
tions of the particles. This was assumed to be proportional to the 
incipient erosion rate. The erosion rate, normalized by the maximum 
rate, is plotted versus the nondimens ional radial distance. For all con- 
figurations and particles the erosion profile is maximum in the area 
.5 <1. The maximum erosion for larger particles occurs near 

— 1, and then drops off sharply with increased radial distance. 

As the particles become smaller a maximum erosion occurs closer to the 
perimeter of the jet nozzle and a second maximum begins to occur at 

occurs after the transition region has begun to 
form and where the onset of turbulence can affect the particles. Fig. 13 
shows the incipient rate of lift entrainment by itself and indicates that, 
in general, the shape of the erosion profile can be attributed to the 
lift forces. 

The results shown in the graphs represent only a few configurations 
selected from the large amount of data collected. In viewing all of the 
results and comparing the many combinations of configurations, it was 
found that the variation of ground erosion rates was principally a func- 
tion of the dimensionless radial distance from the jet centerline and 
the size of the particle for a given nozzle height, (s^ • The ero- 

sion rate was nearly proportional to the particle density except near the 
stagnation point, where the erosion rates were very small. 

From these results it may be concluded that for ground particles of 
less than .04 in. diameter the incipient erosion rate has two maximum 
points. One in the vicinity of the perimeter of the jet nozzle, and the 
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second at a radial distance approximately two diameters from the stagna- 
tion pointo For very small particles and low jet velocities the second 
maximum is the most significant. In all other cases where two maximum 
points exist, the first is the most significant. For ground particles 
whose diameters are greater than .04 in. a single maximum point exists 
near the perimeter of the jet nozzle. As the particle size increases 
the influence of the turbulent region decreases. 

It was concluded that particles of approximately .02 in. diameter 
or less would entrain due to drag even at very high densities or very 
low velocities. For particles larger than approximately .02 in., an 
optimum size of approximately .06 in. exists up to radial distances of 
four nozzle diameters for drag entrainment. An optimum particle size 
most susceptible to entrainment due to lift exists in the range of .07 
to .12 in. diameter for radial distances up to four nozzle diameters. 

Also, it may be concluded that varying the nozzle diameter had very 
little effect on the dimensionless erosion rate; hence in the final analy- 
sis the nozzle diameter was not considered as a parameter. Moreover, 
while the actual erosion rates were approximately linearly dependent upon 
the jet velocity, the dimensionless erosion profile appeared to be inde- 
pendent of jet velocity, with the exception of the order of importance 
of the maximum points. 
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NONDIMENSIONAL RADIAL DISTANCE FROM JET - R/Dn 
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TYPES OP TERRAIN COVERED BY ANALYSIS 
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PLOW FIELD MODEL OP CIRCULAR IMPINGING UNIFORM JET 
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VARIATION OP INVISCID VELOCITY WITH RADIAL STATION 
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VELOCITY PROFILE FOR BOUNDARY LAYER 
NEAR STAGNATION POINT 
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FIGURE U 

GROUND PLANE STATIC PRESSURE DISTRIBUTION - Z/Dn > . 
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VELOCITY PROFILES OF LAMINAR AND TRANSITION 
BOUNDARY LAYER, Z/Dn = .5 
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FIGURE 5b 

VELOCITY PROFILES FOR LAMINAR AND TRANSITION 
BOUNDARY LAYER, Z/Dn = .5 
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VELOCITY PROFILES FOR LAMINAR AND TRANSITION 
BOUNDARY LAYER, Z/Djj-1.0 
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VELOCITY PROFILES OF LAMINAR ANT) TRANSITION 
BOUNDARY LAYER, Z/Djj=1.0 
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VELOCITY PROFILES FOR LAMINAR AND TRANSITION 
BOUNDARY LAYER, Z/Dn=2.0 
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VELOCITY PROFILES OP LAMINAR AND TRANSITION 
BOUNDARY LAYER , Z/Djj = 2.0 
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VELOCITY PROFILES OP LAMINAR AND TRANSITION 
BOUNDARY LAYER , Z/Djj = J 4 . , 0 
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VELOCITY PROFILES FOR LAMINAR AND TRANSITION 
BOUNDARY LAYER , Z/Djj = i;. 0 
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VELOCITY PROFILE OP TURBULENT 
BOUNDARY LAYER - INNER LAYER 
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FIGURE 6b 



VELOCITY PROFILE OF TURBULENT 
BOUNDARY LAYER - OUTER LAYER 
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COMPARISON OF SHEAR LIFT FORCE BY 
CYLINDER AND SPHERE METHODS 
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FIGURE 8a 

CRITICAL PARTICLE SIZE FOR DRAG ENTRAINMENT 
Vj= 65 ft/sec, Z/Dji = l, R/Dn=t.05 
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CRITICAL PARTICLE SIZE FOR DRAG ENTRAINMENT 
Vj 65 ft/sec, Z/Dj| = 1, R/Djj=,1 
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CRITICAL PARTICLE SIZIS FOR DRAG ENTRAINMENT 
Vj= 65 ft/sec, Z/Dii = l, R/Dn-,'2 
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CRITICAL PARTICLE SIZE FOR DRAG ENTRAINMENT 
Vj = 65 ft/sec, Z/Dn=1, R/Dij = .3 
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CRITICAL PARTICLE SIZE FOR DRAG ENTRAINMENT 
Vj = 65 f t/scc , Z/Dii - 1, R/Dn = .1; 
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CRITICAL PARTICLE SIZE FOR DRAG ENTRAIHMENT 
Vj=65ft/3 0C, Z/D,^^-l, R/Dn= .5 
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FIGURE 8g 

CRITICAL PARTICLE SIZE FOR DRAG ENTRAINMENT 
Vj = 65 ft/soc, Z/Dil-1, r/Dn=^2,0 
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CRITICAL PARTICLE SIZE FOR DRAG ENTRAINMENT 
Vj-65 ft/soc, Z/Dn=1, R/Dii=[(.,0 
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CRITICAL PARTICLE SIZE FOR LIFT ENTRAINMENT 
Vj = 6^ rt/soCf Z/Djf‘=^lj R/Dj^ =• ,0^ 



r 




II 

‘ I 



I 



I 

I 

I 

I 

I 

I 

I 

I 



< 



I 




o 

w 

Cd 

o 

M 

Pti 



U9 



CRITICAL PARTICLE SIZE FOR LIFT ENTRAINMENT 
Vj-^ 65 ft/scc, Z/D^^l, R/D^=.1 
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CRITICAL PARTICLE SIZE FOR LIFT ENTRAINMENT 
Vj= 65 ft/sec, Z/Dii=1, R/Di;-.2 
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— i/ft xio-3 
FIGURE 9d 

CRITICAL PARTICLE SIZE FOR LIFT ENTRAINMENT 
Vj==65 ft/scc, Z/Dji^1, R/D]}=:.3 
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CRITICAL PARTICLE SIZE FOR LIFT ENTRAINMENT 
Vj=65 ft/scc, Z/Dn=1, R/Dj>i=.[i. 
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CRITICAL PARTICLE SIZE FOR LIFT EMTRAIMMIINT 
Vj=65 ft/sec, Z/Dii=1, r/D]^=.5 
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FIGURE 9g 

CRITICAL PARTICLE SIZE FOR LIFT ENTRAINMENT 
Vj = 65 ft/sec, Z/Dn = 1, r/Dn = 1.0 
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CRITICAL PARTICLE SIZE FOR LIFT ENTRAINMENT 
Vj = 65 ft/sec, z/Dn=-1, R/Dn-1.33 
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CRITICAL PARTICLE SIZE FOR LIFT ENTRAIIR^IENT 
Vj-65ft/scc, Z/Dn = 1, R/Dn = 2.0 
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CRITICAL PARTICLE SIZE FOR LIFT ENTRAINMENT 
Vj = 65ft/3CC, z/Dn = 1, R/Dn-1i.O 
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CRITICAL PARTICLE SIZE FOR DRAG ENTRAIMENT 
Yj -65 ft/sec, Z/Dk=4,0, R/Di>i = .5 
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CRITICAL PARTICLE SIZE FOR DRx\G ENTRAINI-IEHT 
Vj= 65 ft/sec, Z/Dfi=4.0, R/DJ.J-.667 
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CRITICAL PARTICLE SIZE FOR DRAG ENTRAINMENT 
Vj = 6^ ft/sec, Z/Dn = U.O, R/Dn = .83 
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CRITICAL PARTICLE SIZE FOR DRAG ENTRAINMENT 
Vj=65 ft/sec, Z/Dn= 14 -. 0 , R/Djj = 1.0 
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CRITICAL PARTICLE SIZE FOR DRAG ENTRAINMENT 
Vj=65 ft/sec, Z/Diq=i4..0, R/Djj=1.17 
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CRITICAL PARTICLE SIZE FOR DRAG ENTRAINMENT 
Vj=65 ft/sec, Z/Dj| = U.O, R/Di^=1.33 
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CRITICAL PARTICLE SIZE FOR DRAG ENTRAINMENT 
ft/soc, Z/Djj=4.0, R/Dn-2.0 
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CRITICAL PARTICLE SIZE FOR DRAG ENTRAINMENT 
~ ^5 ft/sec, Z/Djj = ii-.0, R/Dj^=-1|.»0 
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CRITICAL PARTICLE SIZE FOR LIFT ENTRAIlfl^IENT 
Vj -3 65 ft/sec, Z/Djj= 1;.0, R/Dj^xI+.O 
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INCIPIEKT EROSION RATE 
PARTICLE DIAM.- .003 in, Vj= 65 ft/sec, 
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INCIPIENT EROSION RATE 

PARTICLE DIAM.= .007 in, V j = 65 ft/sec, Z/D|^ = 
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INCIPIENT EROSION RATE 
PARTICLE DIAH.= .01 in, Vj- 65 ft/sec, 
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INCIPIENT EROSION RATE 

PARTICLE DIAM.= .02 in, Vj= 65 ft/scc, Z/D 
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INCIPIENT EROSION RATE 
PARTICLE DIAM.= .01^. in, Vj= 65 ft/soc, 
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INCIPIENT EROSION RATE 

PARTICLE DIAM.= .06 in, Vj = 65 ft/sec, Z/Dij = 
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INCIPIENT EROSION RATE 

PARTICLE DIAH, - .125 In, Vj-6^ ft/scc, Z/Djj 
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INCIPIENT EROSION RATE 

PARTICLE DIAM.= .003 in, V j = 500 ft/sec, Z/Djj 
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INCIPIENT EROSION RATE 

PARTICLE DIAM.= .003 in, Vj = $00 ft/sec, Z/Djj 
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INCIPIENT EROSION RATE 
PARTICLE DIAH.= .01 in, V j = 500 ft/soc, 
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INCIPIENT EROSION RATE 

PARTICLE DIAM. = .02 in/ Vj-500 ft/acc, Z/D^ = 
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INCIPIENT EROSION RATE 

PARTICLE DIAM.= .06 in, Vj=500 ft/oec, Z/D^j = 
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INCIPIENT EROSION RATE 

PARTICLE DI.AII, = .125 in, Vj= 500 ft/scc, Z/Dj^ 
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INCIPIENT E:R0SI0N RATE FOR LIFT ENTRAINT-iEI^T 
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INCIPIENT EROSION RATE FOR LIFT ENTRAIWIENT 
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INCIPIENT EROSION RATE FOR LIFT ENTRAINMENT 
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o 







XA 



86 



1.0 1.5 2.0 2.5 3.0 

R/Dn 

FIGURE l[i.f 

INCIPIENT EROSION RATS FOR LIFT ENTRAINMENT 
PARTICLE DIAM.= ,06 in, Vj= 65 ft/occ, ’Z./D}]- 




87 



INCIPIENT EROSION RATE FOR LIFT ENTRAINMENT 
PARTICLE DIAM.= .125 in, Vj=65 ft/scc, Z/Du ^ 
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FIGURE 1 $ 

LAMINAR BOUNDARY LAYER THICKNESS 
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APPENDIX 



FORTRAN PROGRAMS 

LIFT FORCES: 

Input. 

Program is generated from the polynomials describing 
the velocity profiles of the boundary layer, for 
various radial stations* 



MM 

PCOZ 

QSH 

VJ 

DN 

RG 

ZG 

YHVX 


Number of terms of the polynomial 
Coefficients of the polynomial 
Nondimensional shear velocity 
Jet velocity 
Nozzle diameter 

Radial distance from jet centerline 
Nozzle height above ground 

Distance to have maximum velocity in turbulent 
boundary layer 


VMX 


Maximvim velocity in turbulent boundary layer 



Intermediate Results. and Output 



DP 

VI 

YI\ID 

XND 


Particle diameter 
Inviscid velocity 

Nondimensional distance above ground 
Nondimensional roots of velocity profile 
polynomial 


DNS 

DVX 

SDVX 

DUDY 

PLSC 

SLFT 

PLPT 

PLAC 


Particle density 
Derivative of polynomial 
Slope of velocity profile •• 

Slope of velocity profile - 
Shear lift {cylinder method) 

Shear lift (sphere method) 

Lift due to wall constraint 

Acceleration of particle in vortical direction 



Subroutine Root finds the roots of the velocity profile 
polynomi al 
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SUBROUTINE LIFT( MM,PCOZ» OSH , V J , ON , RG , ZG , NCC , YHVX , VMX ) 

DIMENSION DP(10),RP<10)» YGU ( 4 0 > . YGB { 4 0 ) » YL < 70 ) , YG ( 7 0 ) » 

1AR(7o)»VL(7o) »PCOZ(3o ) ,XL(4o),FX(3o)»NLY(5),DVX(30)» 

2 UG6(b),PCOX(10)»R(100),YNU(l20),XND(120)» YNDD(50)»XNDO(50). 

3 VLD(50 ) » YGDD(2o ) 

UP ( 1) = . 06 
0P(2)=. 04 
DP(3)=.q2 
DP(4)=.01 

DP ( 5 j = . 0 07 
DP(6)=. 125 
DP ( 7 ) = . 003 
RGN=RG/DN 
ZGN=ZG/DN 
PRINT 611» DN 

611 FORMAT!/// 22HN0ZZLE DIAMETER (FT) = F5.2 ) 

PRINT 612» VJ 

612 FORMAT! 26HN0ZZLE VELOCITY (FT/SEC) = F7.1 ) 

PRINT 4oo»RGN 

400 F0RMAT!26HN0NDIMEN RADIAL DISTANCE = F20.10) 

PRINT .4 01, ZGN 

401 format ( 28HNONU I men HEIGHTH OF NOZZLE = f20.10) 

DO 115 11=1,7 

PRINT 403, DP! I 1 ) 

4q3 format ! //23HPART I CLE D I AMET ER ( I N . ) = F2o.10) 

RP! 1 I ) = DP! I I)/2. 

GNU=3 . 74bE-7/2. 3769E-3 
IF!NCC-25)3 o2,3o3,3o2 
303 VI=VMX 

GO TO 700 

302 IF!ZG/DN-.5)640,640,66l 
661 KZG = ZG/Ui'I 

GO TO(640, 641,640, 643) »KZG 

640 IF!RGN-.35)60,60,61 
61 IF!RGN-. 75)18,18,19 

19 IF!RGN-1.2)20,20,22 

60 VI=VJ*S0RTF!ABSF!1.88*RGN**2-.225*RGN)) 

GO TO 700 

18 VI=VJ*SQKTF!-2.2*RGN»*2+5.13*rqn-1.26) 

GO TO 700 

20 VI=VJ*SQHTF!-2.68*R6N**2+5.48*RGN-1.82) 

GO TO 700 

22 VI=VJ*SQRTF!2.25*RGN**2-9.45*RGN-^8.96) 

GO TO 700 

641 I F ! RGN- . 35 ) 6q , 6q , 644 

644 IF!RGN-1. 0)645,645,647 

645 VI=VJ*SQRTF!-.D985*RGN**2+1.13*RGN-.241) 

GO TO 700 

647 VI=VJ*SQRTF! . 0 77»rGN** 2- . 64l*RGNt-l . 354 ) 

GO TO 700 

,643 IF!RGN- . 40 )62,62, 654 

654 IF!RGN-1. 2)655, 655, 657 

655 VI=VJ*SQRTF(-.5«RGN**2+1.2*RGN--.3) 
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GO TO 7U0 

657 VI=VJ*SQr<T|- (.0^0‘^*RGl'J**2-.3ia*RGN*.744) 

GO TO 700 

62 VI=VJ*SOHTF ( .625*RGN**2+.o25*RGN) 

700 ANCL=101. 

016 TH=DP( 1 1 ) /ANCL 
NN= ANCL 
NA= ( NR+1 ) /2 
R(1)=WP( 1 I ) 

DO 1 K=2,NA 
AA=K-1 
XLC = AA* TH 

1 R(K)=SQR IF(RP( I I )**2-XLC**2) 

100 IF(NCC-26)901,900,900 
901 IF( NCC-25 > 7 04.705,900 

705 YND( 1)=R(1>/YhVX 

GO TO 104 

704 YNU( 1)=(R( 1)/(12.*DN) )*SQRTF(VJ*DN/GNU) 

GO TO 104 
900 AAA=V1/RG 

YND( l)=R(l>*SORTF(2.*AAA/GNU)/12. 

104 XND(1)=0.0 
lF(NCC-25)508,504,508 

50 4 1 F(YND( 1 )-. 15 )500 , 500,501 
500 MM=8 

DO 511 1Z=1,KM 
511 PCOX( IZ)=PCOZ( IZ) 

GO TO 506 
301 MM=4 

DO 5q2 1=1, mm 

N = 1 8 . 

502 PCOX( 1 )=PCOZ(N) 

506 DO 507 KR=1,MM 
£X1=MM-KW 

FX ( KW ) =PCOX ( KW ) * YND ( 1)’**EX1 

507 XND( 1)=XND( l)-^FX(KW) 

GO TO 673 

508 DO 71 KZ=1,MM 
EX1=MM -KZ+1 

FX(KZ)=PCOZ(KZ)*YND(l)**EXl 
71 XND(1)=XND(1)-"FX(KZ) 

GO TO 673 

105 CALL H00T(1,MM,PC0Z,XND,YND,NCC) 

673 IF(XND(1)-1.0)31,31,32 

32 XND(D = .9999 
31 VL 1 1 ) = V I *X^a)( 1) 

706 ART=0, 

DO 43 JZ=1,MA 
IF( JZ-1)44,44,45 

44 ARE=2.*R( JZ)*3.14159*TH 
GO TO 43 

45 ARE=4.*R( JZ)*3.14159*TH 
43 ART=ART*ARb 

SPA = RP( I I >**2*4. *3. 14159 

ADS=SPA- ART 

ADS=ADS/SPA 

AC0R=1 . +AD5 

FLS = 0 . 

DO 39 IZ=1,NA 
IF(IZ-l) 40,40,41 

40 FL=-TH/2.*2.3769E-3*R(IZ)*VL(1)**2*QSH/144. 
GO TO 39 
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39 FLS = FLS-^FL 

FLSC=FLS*AC0R 

approximate Shear lift of particle by sphere method 

130 YGDD(1)=R<1) 

IF(NCC-25)132»77,129 

77 YNDD(1) = YGUD(1)/YH\/X . . ^ 

GO TO 133 

129 YN0D(1) = YGUD(1)*S0RTF(2.’*AAA/GIMU)/12. 

GO TO 133 

132 YNDD(1)=<YGDD(1)/(12.*un) ) *SQRTF ( V J*DM/GNU ) 

133 SDVX=0.0 

DO 74 KD=1.MM 
MX2=MM-KD 
EX1=MM-KD+1 
I F ( NCC-25 ) 7S , 79, 75 
79 PC02(KD)=PC0X (KD) 

75 IF( KD-MM)180 . 181, 160 

180 DVX(KD)=EX1'»PC0Z(KD)*YNDD(1)**MX2 
GO TO 74 

181 DVX(KD)=PCOZ(KD) 

74 SOVX = SUVX-t-DVX(KD) 

IF(S0VX)182»183,183 

182 IF(NCC-25)184,183,184 
184 SDVX=0.0 

IF(NCC-25)190.183,192 

183 DUDY=SUVx*VMX/YHVX*12 . 

GO TO 193 

192 DUDY=SDVX*S0RTF( 2 . *AAA/GNU) *VI 
GO TO 193 

190 DUDY=SDVX*S0RTF( VJ*DN/GNU) *VI /DN 

193 20N= R(1)/VL<1)*DUL)Y/12, 

CLF = 0 . 345-^0.55 7*ZOiM-^0.8/9*ZON**2 

SLFT=CLF*2.3769E-3*VL(1)**2/2. *3. 14159*R( 1 )**2/144 . 

PRINT 800 

800 format ( 14HNONDIMEN SHEAR , 5X » 1 7HNONO 1 MEN VELOC I T Y , 5X , 17HL I FT DUE TO 
1 SHEaR,5x, 17HLIFT DUE TO SHE AR » , 1 9H A RE A DIFFERENCE 0/0 , 3x , 13 hPR0 
2FILE SLOPE) 

PRINT 801 

801 F0RMaT(3X,8HVEL0CITY,11X,11H0N C YL I NDER , 12X » 1 0 H ( C YL I NDER ) » 13X . 

1 8H(SPHEHE)/) 

PRINT 802,OSH,XND(1) , FLSC , SLF T , ADS , SD VX 

802 format ( 6E20 . lO ) 

GO TO 52 

5 PRINT 200 
200 format ( /5HERROR/ ) 

52 NNT=59 
UN=NNT 

603 dGB ( 1 ) = 0 . 

YH = 0 . 

YLM = NN T-1 
YINC = DP( 1 I )/YLM 
DO 169 J=1»NNT 
IF(J-1)666, 111,112 

111 YG(J)=0. 

GO TO 169 

112 YG( J)=YH+YINC 
YH=YG< J) 

169 CONTINUE 

DO 113 K=2,NNT 
YND( 1 ) =0 . 

IF ( NCC-26 ) 116, 902, 902 



902 YND(K) = YG(»\)*SQRrF{2.*AAA/GNU)/12. 

GO TO 113 

116 IF(NCC-25)300>301,902 
301 YND(K)=YG(K)/YHVX 
GO TO 113 

300 YND(K)=(YG(K)/(12.*0N) ) *SQRTF ( V J*DN/GNU ) 

113 CONTINUE 

671 DO 672 NX=1.NNT 
XND( NX ) = 0 . 

lF(NCC-25)7oa, 720,708 

720 1F(YND(NX)-. 15)721 ,721,701 

721 MM=8 

DO 811 IZ=1,MM 
811 PCOX( 1Z)=PC0Z( IZ) 

GO TO 708 

701 MM =4 

DO 702 1=1, MM 
N = 1 + 8 

702 PCOX{ I )=PCOZ(N) 

DO 707 Kw=l,MM 
EX1=MM-KW 

FX(KW)=PCOX(KW)*YND(NX)**EXl 
• 707 XND(NX)=XND(NX)+FX(KW) 

GO TO 672 
708 DO 37 KL=1,MM 
EX1 = MM-KL-^1 

FX(KL)=PC0Z(KL)*YND(NX)**EX1 

37 XND ( NX ) = XND ( NX ) -^FX ( KL ) 

672 CONTINUE 

DO 95 KC=1,NNT 
IF(XND(KC) >33,35,35 

33 1F(KC-D38,38,51 

38 XND ( KC ) = 0 . 0 
GO TO 35 

5l XND(KC)=XND(KC-1) 

35 1F(XND(KC)-1.0)95,34,34 

34 XND(KC)=.9999 

95 CONTINUE 

FIND AVERAGE VELOCITY 

96 SUMV=0. 

DO 114 L=1,NNT 
VL ( L ) = V I *XND( L ) 

114 SUMV=SUMV+VL(L) 

VLAV=8UMV/UM 

PLIFT=(2.3769E-3/2.)*VLAV**2*3.14167*(3./8,)*RP(II)**2/144, 
PLIFI=-PLIF T 

VOL=4./3.*3.l4l59*KP(II)**3 

DO 730 MO=l,4 

GO TO( 731, 732W33, 734 ) , MO 

731 DNS=75. 

GO TO 735 

732 DNS=90. 

GO TO 735 
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733 dnS=105. 

GO TO 735 

734 DNS=125. 

735 PW I =DNS/1728 . *VOL 
PRINT 736, DNS 

736 FORMAT ( /18HPART I CLE DENSITY = F6.1) 

PLIOT= PL 1 FT-PW r-^FLSC 
PMAS=PWT/32 . 16 

PLAC=PLT0T/PMAS _ 

1 F ( PLA0 842, 843, 843 

842 PLAC=0.0 

843 print 110 

110 format ( /22HAVERAGE LOCAL VE LOC I T Y , 3 X . 1 6HL I F T ON PARTICLE, 3X, 

1 26HNET LIFT FORCE ON PART I CLE , , 21hPaRT I CLE ACCELERATION) 

PRINT 120 

l2o F0HMAT(7X,6H(FT/SEC) ,7X,22 hDUE TO WALL CONSTR A I NT , 30 X , 2lH I N VERTIC 
lAL DIREC [ ION/) 

PRINT 840, VLAV,PLIFT,PLT0T,PLAC 
840 FORMAT ( F20.10,E20.10,5X,F20.10,5X,F20.10/) 

730 CONTINUE 
GO TO 115 

666 PRINT 667 

667 F0KMAT(/5HERR0R) 

GO TO 668 • 

115 CONT INUE 

668 REIURN 
END 

subroutine ROOT(NNi,MM,PCOZ,YND,XND,NCC) 

DIMENSION YND(70),XND(70),PCOZ(30),DvX(30),FX(30) 

DO 210 L=1,NNT 
IF(NCC-24)199, 199,201 
199 IF(YND(L)-. 425 >200,201, 201 
201 XND(L)=0.1 
106 SFX=0. 

SDVX = 0 . 

D0100N=1,MH 

MXl=MM-N+l 

Mx2=mH-N 

EX1=MX1 

FX(N)=PC0Z(N)*XND(L)**MX1 

IF(N-MM)101,102,101 

101 DVX(N)=EX1*PC02(N)*XND(L)*>*MX2 
GO TO 103 

102 DVX(N)=PCOZ(N) 

103 SFX=SFX+FX(N) 

100 SDVX=SDVX+DVX(N) 

RMX=(SFX-YND(L) )/SDVX 
IF(A8SF(RHX)-.01 )108,108,104 

104 XND(L)=XnD(L)-RMX 
MLM=1000/MM 

IF(XND(L)-2.-*MLM)106,204,204 
GO TO 210 . 
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108 IF(NCC-2'^) 109, 109,^^10 

109 IF(MNT-?)^10,210,110 

110 IF(XND(L)-1. 0)210, 203, 203 
203 IF(YNU(L)-2. 0)204, 204, 205 
2o5 XND(L)=.9999 

GO TO 210 

20"^ IF(L-1)^0'^,3o^,3o5 

304 XND(L) =0 . 0 
GO TO 210 

305 XND < L) =XNO I L-1 ) 

GO TO 210 

200 GO TO 210 
210 CONTINUE 
RETURN 
END 
END 



DRAG FORGE: 

Input, 

Program is generated in the same manner as the lift 
program. 

Intermediate Results and Output, 

DELT Boundary layer thickness 

XDEL at top of boundary layer 

FDT Drag due to pressure force 

PPDT Drag due to friciton 

PV7T Particle weight 

PDAC Horizontal acceleration;^ of particle 



SUBROUTINE DRAG(PC0Z,MM,NCC,VJ,RG,ZG,DN,YHVX,VMX) 

DIMENSION DP(10),UN(40),RP(10),YGU(40),YG8(40),YL(70),YG(70), 

1 AR(70),YND(70),XND(70),VL(70),PCOZ(30),XL(4o),FX(30), 

2 XDEL(5),YMAX(5),PCOX(10),ARF(70) 

DPT 1 ) = . 06 

DP(2) = . 04 
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DP(3)=.02 

DP(4)=.o1 

DP(5)=.007 

DP(6)=.i25 

DP(7)=,003 

NN = 3i 

RGN=R6/DN 

ZGN=ZG/DN 

PRINT 611, DN 

611 format ( 22HN0ZZLE DIAMETER (FT) = F5.2) 

PRINT 612, \/J 

612 F0RMAT(26HN0ZZLE VELOCITY (FT/SEC) = F7.1 ) 

PRINT 400/RGN 

400 format ( 26HN0NDIMEN RADIAL DISTANCE = F20.10) 

PRINT 401, ZGN 

401 format (28HN0NDIMEN HEIGHTh OF NOZZLE = F2Q.10) 

DO 115 11=1,7 

PRINT 403, DP( I I ) 

4q 3 F0RMAT( 23hPARTlCLE diameter-in. = F2o,10 > 

RP( I I ) = DP( 1 I )/2. 

UNI =NN 

TK = RP( I I )/(UNI-,5) 

YLC = 0 . 

SUMAR = 0 . 

NA 1 =NN-1 

DO 3 N=1,NAT • 

YL(N)=YLC-^TK/2. 

AM = N 

XL(N)=S0RTK(RP(II)**2-YL(N)**2) 
AR(N)=2.*XL(N)*TK 
ARF ( N) =4 . *XL ( N ) -TK 
I F ( N-1 ) 8 , 5, 8 

5 SUMAR = SUMAR-^3. 1416/2. *ARF(N) 

GO TO 6 

8 SUMAR=SUMAR+3 . 1416-ARF ( N ) 

6 YGU(N)=RP( I I )-^YLC 
YGB(N) = RP( I I )-YLC 

3 YLC=AM*TK 

GNU= . 0000003745/ . 0023769 
IF(NCC-25)302,303,302 
303 VI=VMX 

GO TO 700 

302 IF(ZG/DN-. 5)640, 640, 661 
661 KZG=ZG/DN 

GO TO(640,641,643,643),KZG 

640 1 F ( RGN- . 35 ) 6o , 60 , 61 
61 IF(RGN-.75)16,18,19 

19 1F(RGN-1. 2)20,20,22 

60 VI=VJ*SQRTF(RGN**2+.0175*RGN) 

GO TO 700 

18 VI=VJ*SQRTF(-3.2*RGN**2+5.13*rgN-1.26) 

GO TO 700 

20 VI=VJ*SQRTf(-2.68*RGN**2+5.48*RGN-1.82) 

GO TO 700 

22 Vl=VJ*S0RTF(.256*RGN**2-1.48*RGN-f 2.311) 

GO TO 700 

641 IF ( RGN- . 35 ) 6q , 6o , 644 
64 4 - IF (RGN -1.0 ) 64 5, 6 4 5, 64 7 

645 VI=VJ*SORTF(-. 0985 *RGN**2+1 . 13*RGN- . 24 1 ) 

GO TO 700 

647 VI=VJ*SQRTF(.077*RGN**2-.641*RGN+ 1.354) 

GO TO 700 
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643 IF(RGN-. 40)62. 62,664 

654 IF(RGN-1. 2)655, 655, 657 

655 VI=VJ*SQRTFC-.5*RGN**2+1.2*RGN-.3) 

GO TO 700 

657 V1=VJ*SQRTFC,0^0^*RGN**2-.318*RGN*.744) 

GO TO 7 01) 

62 VI=VJ*SQKTF( . 625.rgn**2+ . o25*rgn ) 

700 iMNT=2*NN-3 
D029N=1 , NNl 
1F(N-NN)24, 26.26 
24 YG(N)=YGB(NN-N) 

GO TO 27 

26 YGCN)=YGU(N-NN+2) 

1F(N-NN)27,27,28 

27 AR(N)=AR(Ni\-N) 

GO TO 977 

28 AR(N)=AR(N-NN+2) 

977 lF(NCC-26)99,900,900 
900 AAA=V1/RG 

YN1)( N) =YG( N) *SQR TF( 2 . *A A A/GNU) /1 2 . 

GO TO 29 

99 IF(NCC-25)978.301,900 
301 YND(N)=YG(N)/YHVX 
GO TO 29 

978 YND(N)=(YG(N)/(12,*DN))*SQRTF(VJ*DN/GNU) 

29 CONTINUE 

672 DO 805 NW=l,NNT 
XND(NW)=0 . O 
lF(NCC-25)608,804,808 

804 IF(YND(NW)-. 15)800, 800, 801 

800 MM=8 

DO 811 IZ=i,MM 
811 PCOX( 1Z)=PC0Z( IZ) 

GO TO 808 

801 MM=4 

DO 802 1=1, MM 

802 PCOX( I )=PCOZ( 1+8) 

806 DO 807 KR=1,MM 
EX1=MM-KW 

FX ( KW ) = PCOX (KW)*YND(NW)**EX1 

807 XND(NW)=XND(NW)+FX(KW) 

GO TO 805 

808 DO 809 KY=1,MM 
EX1=MM-KY+1 

FX(KY ) = PCOZ(KY)*YN'D(NW)**EXl 

809 XND(NR)=XND(NW)+FX(KY) 

805 CONTINUE 

75 DO 32 KC=1,NNT 
1F(XND(KC) )33,670,670 

33 IF(KC-1)72,72,73 

72 XND(KC)=0.0 
GO TO 670 

73 XND(KC)=XND(KC>1) 

670 IF ( XND ( KC ) -1 . 0 ) 32, 32,34 

34 XND(KC)=.9999 
32 CONTINUE 

1F( I 1 -1) 74,74,80 

74 IF(NCC-25)78,76, 79 

79 DELT=4, 6*12. /SQRTF(2.*AAA/ GNU) 

GO TO 80 

76 XDEL=0.99 

CALL ROOTd, MM, PCOX, XDEL, YMAX,NCC) 
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UFLT=YMAX ( 1 ) * YHVX 
GO TO 80 
78 XDbL(l)=.99 

CALL RUOl (lWlM,PCOX.XDEL.YhAX,NCC) 

DELT = YMAX(1 ) »UN/S0R IT(VJ*DN/GNU)*12. 

60 GO TO 31 
31 D035K=l,NNi 

35 VL(K)=V1*XiMO(K) 

82 NMID=(NM+l)/2 

RN=VL(NMJU)*'DP( I I )/(l2, *GNU) 

IF(RN-100. )730,730,731 

731 1F(RN-10000. )732,732,733 
730 CD=30 . / ( RN** . 699 ) 

GO TO 734 

732 CD=1 .44/(RN**.o396) 

GO TO 734 

733 CD=l.l8 

734 FOr=0. 

D036K = 1» NNT 

rD=2.3769E-3*CD*VL(K)**2*AR(K)/288. 

36 FDI=FDT+FD 
ARS=4.*3.14159*RP(II)**2 
ARD=ARS-8UhAR 

ARD = ARlj/ARS 
ACOR = l , -^ARU 
CF= . 664/SORTF (RN ) 

FFDT=0. 

DO40M=l, NNl 

FFD=2 . 3769b-3*CF*VL (M) **2*ARF ( M) /288 . 

40 FFDT=FFD r+FFD 
CF8=4 . *CF 

rOL=ABSF(CFB*. 00001 ) 

CALL CUBE (CFB» TOL, CFB3 ) 

CBU=0 . 1/CFB3 
FBDT=0 . 

UO 42 MD=1,NNT 

FBD=2.3769E-3*CB0*VL(MU)**2*AR(MU)/288. 

42 F6UT=FBDT+FBD 
PRINT 38 

38 F0RMAT(/4X, 16HPRESSURE DRAG-LB, 4X , 16HFR 1 CT I ON DR A G - L B , 8X , 1 2HB A SE I 
1RAG-LB,2X,24hB 0UNDARY LAYER TH I CKNESS , 2X , 3l HSPHERE APPROXIMATION I 
21FFERENCE// ) 

PRINT 41, FUT,FFDT,FBDT,DELT, ARD 

41 F0RMAT( 3E20.10,2F20.10) 

VOL=4./3.*3.l4l59*RP(II)**3 
DO 720 M0=l,4 

GO T0(721,722,723,736),M0 

721 DNS=75. 

GO TO 724 

722 DN8=90. 

60 TO 724 

723 DNS=105. 

GO TO 724 

736 DNS=125. 

724 Pw1=DNS/1728. *VOL 
PRINT 725, DNS 

725 F0RMAT(/18HPARTIClE DENSITY = F6.D 
TOIL=(FD1+FFDT-,707*PWT)*ACOR 
PRINT 620, TOTL 

620 format ( /34HNET FORCE IN DIRECTION OF MOTION = E20.10) 

PMAS=PWT/32 . 16 
PDAC=TOTl/PMaS 
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IF(PDAC)a42/843,843 
842 PDAC=0.0 
8^3 PRINT 621, PDAC 

621 F0HMAT(/56HPARTICLfc ACCELERATION IN HORIZONTAL DIRECTION (FT/SEC) 
1= F20.10) 

720 CONTINUE 
115 CONTINUE 

IF(NCC-24)668.668,599 
599 CONTINUE 
668 REIURN 
END 
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